{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting Temporal changes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will plot all heights regardless of locations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One track would be plotted as an example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import h5py\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "import datetime\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "import statistics\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/download/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.h5'))\n",
    "## choose a file as an example\n",
    "\n",
    "##### function of reading alt13\n",
    "def alt13_to_df(filename, beam):    \n",
    "    f = h5py.File(filename, 'r')\n",
    "    f_beam = f[beam]\n",
    "    lat = f_beam['segment_lat'][:]\n",
    "    long = f_beam['segment_lon'][:]\n",
    "    ws = f_beam['ht_water_surf'][:]\n",
    "    ws_sd = f_beam['stdev_water_surf'][:]\n",
    "    ws_err = f_beam['err_ht_water_surf'][:]\n",
    "    ortho = f_beam['ht_ortho'][:]\n",
    "    wd = f_beam['water_depth'][:]\n",
    "    alt13_df = pd.DataFrame({'Latitude':lat,'Longitude':long,'SurfaceH':ws,\n",
    "                            'SH_SD':ws_sd, 'SH_error':ws_err,'OrthoH':ortho,\n",
    "                            'WaterD':wd})\n",
    "    return alt13_df\n",
    "\n",
    "### We set 'gt2l' as a beam of example for tracking\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for ff in files:\n",
    "    try:\n",
    "        D_dict[ff]=alt13_to_df(ff, 'gt2l')\n",
    "    except KeyError as e:\n",
    "        ##print(f'file {ff} encountered error {e}')\n",
    "        error_count += 1\n",
    "        \n",
    "### bounds of Tonle Sap Lake\n",
    "sp_ex = [103.643, 104.667, 12.375, 13.287]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot of one track"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### I need to make a col and draw it repetitively"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for record the beam, date and RGT, as well as other information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def alt13_to_df_beams(filename):    \n",
    "    f = h5py.File(filename, 'r')\n",
    "    rgt = str(filename).split(\"_\")[5][0:4]\n",
    "    cycle = int(str(filename).split(\"_\")[5][4:6])\n",
    "    version = int(str(filename).split(\"_\")[5][6:8])\n",
    "    ymd = str(filename).split(\"_\")[4][0:8]\n",
    "    ymd_trans = datetime.datetime(int(ymd[0:4]),int(ymd[4:6]),int(ymd[6:8]))\n",
    "    date = ymd_trans.strftime(\"%Y-%m-%d\")\n",
    "    beam_lst = list(f)[2:-1]\n",
    "    alt13_df = pd.DataFrame()\n",
    "    for beam in beam_lst:\n",
    "        f_beam = f[beam]\n",
    "        lat = f_beam['segment_lat'][:]\n",
    "        long = f_beam['segment_lon'][:]\n",
    "        ws = f_beam['ht_water_surf'][:]\n",
    "        ws_sd = f_beam['stdev_water_surf'][:]\n",
    "        ws_err = f_beam['err_ht_water_surf'][:]\n",
    "        ortho = f_beam['ht_ortho'][:]\n",
    "        wd = f_beam['water_depth'][:]\n",
    "        df_beam = pd.DataFrame({'Beam': beam ,'RGT':rgt,'Cycle': cycle, 'Date':date, 'Date_num':int(ymd),'Ver.':version,\n",
    "                                 'Latitude':lat,'Longitude':long,'SurfaceH':ws,\n",
    "                                'SH_SD':ws_sd, 'SH_error':ws_err,'OrthoH':ortho,\n",
    "                                'WaterD':wd})\n",
    "        alt13_df = alt13_df.append(df_beam, ignore_index = True)\n",
    "        \n",
    "    return alt13_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_beams = alt13_to_df_beams(files[37])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5508cf4abc994b378c20b63c8820c448",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(test_beams['Date'],test_beams['SurfaceH'],'.',markersize=0.25, label='all segements')\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot of all tracks with all beams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "a_tracks = pd.DataFrame()\n",
    "for ff in files:\n",
    "    a_tracks = a_tracks.append(alt13_to_df_beams(ff), ignore_index = True)\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Sorting to draw graph\n",
    "a_tracks.sort_values(by=['Date_num'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "312fad71f4634e0f9625e19e797715f5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "ax.plot(a_tracks['Date'],a_tracks['SurfaceH'],'.',markersize=0.25, label='all segements')\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the tracks by grouping their locations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Map all tracks with labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drawing a just track first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "75a97a59e2524c60a0d5d767cc084a68",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Tracks on TSL')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.set_extent(sp_ex, crs=ccrs.PlateCarree())\n",
    "plt.scatter(a_tracks['Longitude'], a_tracks['Latitude'],s=1)\n",
    "\n",
    "request = cimgt.Stamen('terrain-background')\n",
    "ax.add_image(request, 10)\n",
    "plt.title(\"Tracks on TSL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot water levels according to cycles.\n",
    "\n",
    "Before it, I will plot the third track, and see the temporal gaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Select the 3rd track\n",
    "track_3rd = a_tracks.loc[a_tracks['Cycle'] == 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "98e56f6700da4e9b8ef4f6fffed860f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "NameError",
     "evalue": "name 'track_3rd' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-a38288607efe>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprojection\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mccrs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPlateCarree\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_extent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msp_ex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcrs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mccrs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPlateCarree\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrack_3rd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Longitude'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrack_3rd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Latitude'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mrequest\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcimgt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mStamen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'terrain-background'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'track_3rd' is not defined"
     ]
    }
   ],
   "source": [
    "#### Mapping\n",
    "fig = plt.figure(figsize=(6,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.set_extent(sp_ex, crs=ccrs.PlateCarree())\n",
    "plt.scatter(track_3rd['Longitude'], track_3rd['Latitude'],s=1)\n",
    "\n",
    "request = cimgt.Stamen('terrain-background')\n",
    "ax.add_image(request, 10)\n",
    "plt.title(\"Tracks on TSL\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(track_3rd['Date_num'].max()-track_3rd['Date_num'].min())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is about 201 day. Too long to see the water changes.\n",
    "\n",
    "I will see some papers and see how they make a graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are my plans to do today.\n",
    "1. Plot outliers\n",
    "2. Add averaged points of the day\n",
    "3. See the relationship of water level with extent\n",
    "4. Newly apply to Bugon Res."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot outliers and averaged heights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outliers\n",
    "IQR = q0.75 − q0.25\n",
    "LSHoutlier > q0.75 + 1.5 ∗ IQR or LSHoutlier < q0.25 − 1.5 ∗ IQR\n",
    "\n",
    "I need to do it by beams, so I need to modify the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_tracks.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function for removing outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "### outliers\n",
    "def out_rmv(df,var):\n",
    "    ### lower (25%)\n",
    "    lwq = df[var].quantile(0.25)\n",
    "    ### upper (75%)\n",
    "    hwq = df[var].quantile(0.75)\n",
    "\n",
    "    ### IQR\n",
    "    iqr = hwq - lwq\n",
    "\n",
    "    lw_out = lwq-1.5*iqr\n",
    "    hw_out = hwq+1.5*iqr\n",
    "\n",
    "    ###LSHoutlier > q0.75 + 1.5 ∗ IQR or LSHoutlier < q0.25 − 1.5 ∗ IQR\n",
    "    return df.loc[(df[var] >= lw_out) & (df[var] <= hw_out)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function for stacking DB without outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def alt13_to_df_beams_out(filename):    \n",
    "    f = h5py.File(filename, 'r')\n",
    "    rgt = str(filename).split(\"_\")[5][0:4]\n",
    "    cycle = int(str(filename).split(\"_\")[5][4:6])\n",
    "    version = int(str(filename).split(\"_\")[5][6:8])\n",
    "    ymd = str(filename).split(\"_\")[4][0:8]\n",
    "    ymd_trans = datetime.datetime(int(ymd[0:4]),int(ymd[4:6]),int(ymd[6:8]))\n",
    "    date = ymd_trans.strftime(\"%Y-%m-%d\")\n",
    "    beam_lst = list(f)[2:-1]\n",
    "    alt13_df = pd.DataFrame()\n",
    "    for beam in beam_lst:\n",
    "        f_beam = f[beam]\n",
    "        lat = f_beam['segment_lat'][:]\n",
    "        long = f_beam['segment_lon'][:]\n",
    "        ws = f_beam['ht_water_surf'][:]\n",
    "        ws_sd = f_beam['stdev_water_surf'][:]\n",
    "        ws_err = f_beam['err_ht_water_surf'][:]\n",
    "        ortho = f_beam['ht_ortho'][:]\n",
    "        wd = f_beam['water_depth'][:]\n",
    "        df_beam = pd.DataFrame({'Beam': beam ,'RGT':rgt,'Cycle': cycle, 'Date':date, 'Date_num':int(ymd),'Ver.':version,\n",
    "                                 'Latitude':lat,'Longitude':long,'SurfaceH':ws,\n",
    "                                'SH_SD':ws_sd, 'SH_error':ws_err,'OrthoH':ortho,\n",
    "                                'WaterD':wd})\n",
    "        alt13_df = alt13_df.append(df_beam, ignore_index = True)\n",
    "    \n",
    "    alt13_df_out = pd.DataFrame()\n",
    "    for beam_type in beam_lst:\n",
    "        df_select = alt13_df.loc[alt13_df['Beam']==beam_type]\n",
    "        alt13_df_out = alt13_df_out.append(out_rmv(df_select,'SurfaceH'), ignore_index = True)\n",
    "        \n",
    "    return alt13_df_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def alt13_a_out(db,var2):\n",
    "    alt13_df = pd.DataFrame()\n",
    "    date_list = db['Date_num'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date_num']==prd]\n",
    "        beam_lst = db['Beam'].unique()\n",
    "        for beam_type in beam_lst:\n",
    "            df_select = s_db.loc[s_db['Beam']==beam_type]\n",
    "            alt13_df = alt13_df.append(out_rmv(df_select,var2), ignore_index = True)\n",
    "    return alt13_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_tracks_out=alt13_a_out(a_tracks,'SurfaceH')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared the numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(a_tracks_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(a_tracks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function for making the averaged heights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(statistics.mean(test['SurfaceH']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(test['Date'].unique()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for averaginf the height\n",
    "def alt13_mean(db,var2):\n",
    "    mean_df_list = pd.DataFrame()\n",
    "    date_list = db['Date_num'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date_num']==prd]\n",
    "        av_height = statistics.mean(s_db[var2])\n",
    "        mean_df = pd.DataFrame({'Av_level':av_height,'Date':s_db['Date'].unique(), 'Date_num':s_db['Date_num'].unique()})\n",
    "        mean_df_list = mean_df_list.append(mean_df, ignore_index = True)\n",
    "    return mean_df_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "av_h = alt13_mean(a_tracks_out,'SurfaceH')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "15e16369f9684bcc8a5b8ee682220cb5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "ax.plot(av_h['Date'],av_h['Av_level'],'+',markersize=4, label='all segements')\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot segments and segments without outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4549db9687cd4855b74b16f50e0ab073",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(10,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "ax.set_ylim([-20,0])\n",
    "ax.plot(a_tracks['Date'],a_tracks['SurfaceH'],'.', color = 'orange', markersize=0.25, label='all segements')\n",
    "ax.plot(a_tracks_out['Date'],a_tracks_out['SurfaceH'],'.', color = 'green',markersize=0.4, label='Without outliers')\n",
    "ax.plot(av_h['Date'],av_h['Av_level'],'+', color = 'red',markersize=3, label='Averaged')\n",
    "ax.plot(av_h['Date'],av_h['Av_level'],color = 'red',linewidth=0.25)\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['2018-10-28', '2018-11-01', '2018-11-02', '2018-11-30',\n",
       "       '2018-12-01', '2018-12-30', '2019-01-03', '2019-01-27',\n",
       "       '2019-01-31', '2019-02-01', '2019-03-01', '2019-03-30',\n",
       "       '2019-04-03', '2019-05-02', '2019-05-31', '2019-07-28',\n",
       "       '2019-08-01', '2019-08-29', '2019-08-30', '2019-09-28',\n",
       "       '2019-10-27', '2019-10-31', '2019-11-28', '2019-11-29',\n",
       "       '2019-12-27', '2019-12-28', '2020-01-01', '2020-01-25',\n",
       "       '2020-01-29', '2020-01-30', '2020-02-23', '2020-02-27',\n",
       "       '2020-02-28', '2020-03-27', '2020-03-28'], dtype=object)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_tracks['Date'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load changes of water extents of the TSL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I used GEE to extract the water level changes using MODIS. I made monthly water mask for the cloud-free images and extracted the extents. For filling some gaps due to the clouds, images of one and two years before are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/jovyan/ICESat_water_level/extraction'"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "we = pd.read_csv(\"./data/200706_water_TSL.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Create the col of readable date\n",
    "date_lst = []\n",
    "for dates in we['Date']:\n",
    "    ymd_trans = datetime.datetime(int(str(dates)[0:4]),int(str(dates)[4:6]),int(str(dates)[6:8]))\n",
    "    date_trans = ymd_trans.strftime(\"%Y-%m-%d\")\n",
    "    date_lst.append(date_trans)\n",
    "we.insert(2, \"Date2\", date_lst, True) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    system:index         Area       Date2  id      Date\n",
      "0              0  3370.644055  2018-09-01   0  20180901\n",
      "1              1  3703.477608  2018-10-01   1  20181001\n",
      "2              2  2012.177456  2018-11-01   2  20181101\n",
      "3              3  1836.499641  2018-12-01   3  20181201\n",
      "4              4  2288.920053  2019-01-01   4  20190101\n",
      "5              5  2287.976726  2019-02-01   5  20190201\n",
      "6              6  2124.563827  2019-03-01   6  20190301\n",
      "7              7  2214.760226  2019-04-01   7  20190401\n",
      "8              8  2559.180276  2019-05-01   8  20190501\n",
      "9              9  2575.179994  2019-06-01   9  20190601\n",
      "10            10  2294.466438  2019-07-01  10  20190701\n",
      "12            12  3291.718026  2019-09-01  12  20190901\n",
      "13            13  2508.601048  2019-10-01  13  20191001\n",
      "14            14  2715.205733  2019-11-01  14  20191101\n",
      "15            15  2334.157314  2019-12-01  15  20191201\n",
      "16            16  2348.783952  2020-01-01  16  20200101\n",
      "17            17  2869.605395  2020-02-01  17  20200201\n",
      "18            18  2904.915253  2020-03-01  18  20200301\n",
      "19            19  2401.990595  2020-04-01  19  20200401\n"
     ]
    }
   ],
   "source": [
    "## remove 190801 due to spurring values\n",
    "we = we.drop(11)\n",
    "print(we)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "341d744a247e4f3da5c043702996a681",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### plot\n",
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "\n",
    "x = we['Date2']\n",
    "y = we['Area']\n",
    "plt.plot(x,y, linewidth=1,color = '#45b6fe')\n",
    "ax.plot(x,y,'o',markersize=4, label='MODIS-Driven Extents',color = '#1c4966')\n",
    "\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Extents')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Extents, km^2')\n",
    "fig.autofmt_xdate()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Process the data for drawing the graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing Join in python is difficult for me so I processed it via R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "##3 Save heights\n",
    "av_h.to_csv(\"./data/200706_waterlevel_TSL.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "joined = pd.read_csv(\"./data/joined2.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Unnamed: 0      Date2   Av_level         Area\n",
      "0             1   9/1/2018        NaN  3370.644055\n",
      "1             2   9/2/2018        NaN          NaN\n",
      "2             3   9/3/2018        NaN          NaN\n",
      "3             4   9/4/2018        NaN          NaN\n",
      "4             5   9/5/2018        NaN          NaN\n",
      "..          ...        ...        ...          ...\n",
      "574         575  3/28/2020 -14.823683          NaN\n",
      "575         576  3/29/2020        NaN          NaN\n",
      "576         577  3/30/2020        NaN          NaN\n",
      "577         578  3/31/2020        NaN          NaN\n",
      "578         579   4/1/2020        NaN          NaN\n",
      "\n",
      "[579 rows x 4 columns]\n"
     ]
    }
   ],
   "source": [
    "print(joined)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the gaps, the graph is not drawn well. I will interpolate them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0      3370.644055\n",
      "1              NaN\n",
      "2              NaN\n",
      "3              NaN\n",
      "4              NaN\n",
      "          ...     \n",
      "574            NaN\n",
      "575            NaN\n",
      "576            NaN\n",
      "577            NaN\n",
      "578            NaN\n",
      "Name: Area, Length: 579, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(joined['Area'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a09dd52304274730bb01d486bba980f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### plot\n",
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "\n",
    "x = joined['Date2']\n",
    "y = joined['Area']\n",
    "ax.plot(x,y,'o',markersize=4, label='MODIS-Driven Extents',color = '#1c4966')\n",
    "\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Extents')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Extents, km^2')\n",
    "fig.autofmt_xdate()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d2182eea8954431870d248712a3854f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "x= joined['Date2']\n",
    "y1= joined['Area']\n",
    "y2= joined['Av_level']\n",
    "\n",
    "color = 'tab:red'\n",
    "ax1.grid()\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('Water Extents', color=color)\n",
    "ax1.plot(x,y1,'x',markersize=2, label='MODIS-Driven Extents',color = color)\n",
    "##ax1.plot(x,y1,'o',markersize=4, label='MODIS-Driven Extents',color = color)\n",
    "ax1.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax2.set_ylabel('Water Level', color=color)  # we already handled the x-label with ax1\n",
    "##ax2.set_ylim([-20,0])\n",
    "ax2.plot(x,y2,'+',markersize=2, label='MODIS-Driven Extents',color = color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "fig.tight_layout()  # otherwise the right y-label is slightly clipped\n",
    "fig.autofmt_xdate()\n",
    "ax1.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install GEE to process the water extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does not work well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install earthengine-api"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install earthengine-api --upgrade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://accounts.google.com/o/oauth2/auth?client_id=517222506229-vsmmajv00ul0bs7p89v5m89qs8eb9359.apps.googleusercontent.com&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fearthengine+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fdevstorage.full_control&redirect_uri=urn%3Aietf%3Awg%3Aoauth%3A2.0%3Aoob&response_type=code&code_challenge=78G-PZM_TXbWeyeicjy7XoU9QMf-zRcpMH4q5CNt7Qg&code_challenge_method=S256>https://accounts.google.com/o/oauth2/auth?client_id=517222506229-vsmmajv00ul0bs7p89v5m89qs8eb9359.apps.googleusercontent.com&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fearthengine+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fdevstorage.full_control&redirect_uri=urn%3Aietf%3Awg%3Aoauth%3A2.0%3Aoob&response_type=code&code_challenge=78G-PZM_TXbWeyeicjy7XoU9QMf-zRcpMH4q5CNt7Qg&code_challenge_method=S256</a></p>\n",
       "        <p>The authorization workflow will generate a code, which you\n",
       "        should paste in the box below</p>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Enter verification code:  4/1gEBxgk5Mq8dmKFQrrkfsFNF4ilSoLvEQ06SuhE2M0u9rmwFn7paSps\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Successfully saved authorization token.\n"
     ]
    }
   ],
   "source": [
    "ee.Authenticate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi = ee.ImageCollection(\"LANDSAT/LC08/C01/T1_8DAY_NDWI\").filterDate('2018-01-01','2020-07-01')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'properties'"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(ndwi.getInfo())[4]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
